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ABSTRACT 

We  discuss  techniques  for  the  estimation  of  nonlinearities  and  state-dependent 
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INTRODUCTION 

In  this  presentation  we  discuss  methods  for  inverse  or  parameter  estimation  problems 
which  can  be  employed  as  quantitative  modeling  techniques  in  models  for  distributed 
(spatially,  age,  size,  etc.)  biological  systems.  In  this  context  they  may  be  useful  in  attempts 
to  understand,  elaborate  on,  or  further  refine  details  of  specific  mechanisms  for  dispersal, 
growth,  interaction,  etc.  in  wide  classes  of  models.  We  have  also  used  these  techniques  in  a 
number  of  biologically  related  problems  [1]  such  as  bioturbation  [12],  [14],  [15]  and 
climatology  [19].  In  addition  to  an  overview  of  ideas  underlying  these  techniques,  we  shall 
present  here  brief  discussions  and  some  findings  on  two  specific  biological  problems  for 
which  we  are  currently  using  them  successfully. 

A  typical  inverse  problem  entails  some  given  or  hypothesized  dynamical  model  with 
"parameters"  q  (often  temporally  and/or  "spatially"  or  even  state  dependent)  and  "states" 
u(t.x;q),  0$t<T,  xcO,  which  depend  on  the  parameters  through  a  dynamical  system  of 

A 

equations.  One  has  observations  or  data  u^j  for  u(tj,Xj;q)  and  wishes  to  choose,  from  some 
admissible  parameter  set  Q,  parameters  q  so  as  to  give  a  best  fit  of  the  model  to  the  data. 
For  example,  we  might  have  a  hypothesized  model  for  transport 

a*  k  I;)' 

with  initial  and  boundary  conditions  also  possibly  depending  on  the  unknown  parameters 
q=(l'.f>.F).  Given  data  u^j,  we  seek  to  minimize  a  fit>to-data  criterion  such  as  a  least-squares 

J(q)  =  E  |u(tj,Xj;q)  -  Ujjl*  (2) 

■J 

over  a  specified  class  Q  of  functions  (I'.D.F)  so  as  to  obtain  a  best  estimate  q  =  (P.D.F).  In 
addition  to  obtaining  estimates  for  q,  usually  one  desires  to  analyze  in  some  way  the  "goodness" 
of  the  model  in  describing  the  phenomena  one  is  modeling.  We  shall  elaborate  on  some  related 
questions  in  this  regard  below. 


The  methods  we  discuss  briefly  here  can  be  powerful  modeling  tools  when  carefully 
and  correctly  used.  Some  of  the  novel  features  of  our  recent  efforts  include  the 
capabilities  for  estimation  of  (i)  state  or  density  dependent  dispersal  coefficients  such  as  D 
above  in  (I),  (ii)  system  nonlincarities  such  as  F  in  (1),  and  (iii)  boundary  parameters  in 
both  simple  and  not-so-simpic  boundary  conditions  (we  give  an  example  in  the  discussions 
on  size  dependent  models  below).  Furthermore,  there  arc  a  number  of  modeling  related 
questions  that  one  might  hope  to  address  from  a  theoretical  or  computational  (or  both) 
viewpoint  with  the  aid  of  these  techniques.  These  include: 

(a)  Experimental  design  (3J,  (4J.  [5]:  What  is  the  appropriate  data  required  to  support  analysis 
of  a  particular  model  or  mechanism?  E.g.,  How  many  time  vs.  spatial  observations  must  be 
made,  or  what  type  of  initial  data  is  needed  to  study  movement  patterns? 

(b)  Robustness  of  model  parameters  [I],  (3],  (SJ,  (I6J,  (I9J:  Do  the  problem  formulation  and 

the  methods  enjoy  certain  stability  properties?  E.g.,  Do  the  parameter  estimates  and  the 
estimation  methods  depend  continuously  on  observation  noise,  initial  data,  amounts  of 
data  available,  problem  constraints,  etc.? 

(c)  Idcntifiabilitv  [18],  (25]:  Is  the  map  from  the  parameter  space  to  the  observation 

sufficiently  well-behaved  so  that  the  methods  can  produce  unique  estimates? 

(d)  Model  comparison  [3],  (5),  [19]:  Can  one  make  evaluations  regarding  the  importance 

and  type  of  mechanisms  needed  to  model  given  phenomena?  E.g.,  Which  is  more 
important  in  particular  transport  phenomena:  convection,  diffusion,  nonlinear  effects, 
dynamic  (time  varying)  vs.  heterogeneous  (spatially  varying)  terms  or  coefficients?  What 
level  of  refinement  in  modeling  terms  can  be  supported  by  the  experimental  design  and 
data?  Do  model  refinements  yield  statistically  significant  improvements  in  explanation  of 
the  data? 

CONCEPTUAL  CONSIDERATIONS 

Wc  next  outline  briefly  certain  ideas  related  to  the  problems  and  methods  that  arc  the  focus 
of  this  presentation.  At  the  same  time  wc  shall  indicate  some  questions  that  may  arise  in  either 
theoretical  or  computational  aspects  of  investigations  using  the  methods.  These  discussions  can 
be  made  precise  and  mathematically  rigorous,  but  for  the  sake  of  brevity,  wc  shall  not  do  that 
here. 

For  the  purposes  of  illustration,  we  return  to  the  problem  of  minimizing  the  functional  of 
(2)  subject  to  the  system  (1)  relating  the  states  and  parameters.  Such  problems  lead  to  the  need 
for  optimization  techniques  for  constrained  problems  that  are  infinite  dimensional  in  nature. 
One  has  a  system  with  states  (t,x)-'u(t,x)  in  some  infinite  dimensional  function  space  X  and 
parameters  (t,x)-^(t,x),  or  (t,x,u)-^i(t,x,u)  if  the  parameters  arc  state  dependent,  in  some  infinite 
dimensional  function  space  Q.  These  problems  can  be  concisely  stated  in  a  theoretical 
framework  using  either  the  theory  of  semigroups  or  evolution  operators,  or  the  theory  of 
scsquilincar  forms  in  Hilbert  spaces.  Wc  won’t  pursue  the  details  here,  but  refer  the  reader  to 
[21,  [7],  (9],  (13],  (15]. 

In  any  case,  this  leads  to  the  recognition  that  in  order  to  effectively  develop 
computational  techniques,  one  must  introduce  approximation  schemes  for  the  state  and 


parameter  spaces.  That  is,  one  needs  families  of  finite  dimensional  spaces  and  (such 
as  finite  elements,  splines,  spectral  families)  such  that  approximates  X  well  as  N-**  and 

approximates  Q  as  M-".  (We  shall  (imprecisely)  write  this  as  X^-*X,  or  simply 

N-*“  M-*“,  in  the  discussions  below.)  One  then  must  develop  schemes  to  solve  the 

approximate  problems  obtained  when  u  in  J  of  (2)  is  replaced  by  the  approximate  states 
u’^tX'^  satisfying  some  equation  approximating  (1).  Minimization  is  carried  out  over 
yielding  approximate  best-fit  parameters  Thus,  the  algorithms  we  have  developed 

and  used  (e.g.  see  [3],  [5],  [8]  for  details)  entail  iterative  optimization  techniques  combined 
with  appropriately  chosen  approximation  schemes  based  on  families  X*^, 

Among  the  important  questions  associated  with  these  approximation  ideas  are  those  of 
method  convergence  and  method  stability.  In  the  first,  one  must  argue  that  as  X^-X, 

Q^-O,  where  q  is  a  solution  to  the  original  problem  involving  (1)  and  (2).  That  is,  one  must 

assure  fidelity  of  the  estimates  under  sufficiently  accurate  approximation  of  state  and 
parameter  spaces.  The  concept  of  stability  is  related  to  a  continuous  dependence  of  the 
estimates  on  the  observed  data,  X*^,  and  QH  More  precisely,  if  q^’^(u)  denotes  solutions  to  the 
approximate  problems  corresponding  to  state  space  X^,  parameter  space  and  data  u,  and  if 
(u*^)  is  a  sequence  of  data  with  u*^-*u,  then  one  desires  to  guarantee  that  q^’  *^(u*^)-^(u)  as 

_  ^  A 

N,M,K-**,  where  q(u)  is  a  solution  of  the  original  problem  with  data  u.  That  is,  fidelity  of  the 
estimates  will  hold  with  sufficiently  small  noise  in  the  observations  as  well  as  sufficiently 
accurate  approximations  of  the  state  and  parameter  spaces.  For  further  discussions  see  [1], 
[18]. 

One  can  develop  a  general  theory  to  guarantee  theoretically  and  computationally 
well-behaved  algorithms  based  on  the  ideas  we  have  used  in  a  wide  class  of  problems  including 
the  examples  discussed  below.  The  arguments  rely  heavily  on  ideas  from  functional  analysis, 
approximation  theory  and  compactness.  We  refer  the  reader  to  (1),  (4),  [8],  [15],  [18],  [25]  for 
further  elaboration  and  details.  We  only  note  here  that  fundamental  to  all  these  convergence 
and  stability  results  is  the  establishment  that  u*^(t,x;q^)-u(t,x;q)  in  some  sense  (i.c.,  in  an 
appropriate  X-topology)  whenever  q^-^i  in  an  appropriate  sense  (i.e.,  in  a  Q-topology).  For 
further  discussion  of  mathematical  ideas,  and  implementation  and  testing  of  the  methods,  we 
refer  to  the  presentations  in  [3],  [5]  in  addition  to  those  references  cited  above.  Here  we 
discuss  several  projects  in  which  these  methods  are  playing  a  fundamental  role  and  outline  sonic 
new  results  in  two  areas. 


We  have,  in  collaboration  with  P.  Kareiva  (U.  Washington),  considered  a  number  of  aspects 
of  insect  movement  and  growth.  In  several  cases  our  quantitative  methods  have  proved  useful 
in  planning  the  experiments  as  well  as  actually  investigating  various  models.  Among  the 
investigations  we  have  pursued  are: 

(i)  quantifying  "initial  disturbance”  effects  in  dispersal  rates  for  flea  beetle  movement  in 
mark-recapture  experiments  in  cultivated  collard  patches  [3],  [5],  [21],  [22],  [23]; 

(ii)  studying  the  effects  of  density-dependent  dispersal  rates,  nonlinear  growth,  interaction, 
and  predation  in  multiple  species  models  such  as  those  for  ladybug-aphid-goldenrod 


experiments  [10],  [24]; 

(iii)  quantifying  'preferred  direction*  components  in  cabbage  root  fly  movement  in 
two-dimensional  domains  [6],  [31]. 

In  regard  to  the  investigations  of  (ii),  we  note  that  the  methods  can  be  used  effectively  to 
estimate  the  shape  of  density-dependent  dispersal  coefficients  D  and  nonlinear  growth  terms  f 
in  models  of  the  form 


where  D(t,x,  ■)  has  the  form  depicted  in  Figure  1. 


Figure  1 

Such  dispersal  coefficients  represent  a  rate  that  is  bounded  below  and  above  (basal  and 
saturation  limits)  and  depends  linearly  on  the  density  between  these  bounds.  Problems  with 
density-dependent  dispersal  have  received  attention  elsewhere  [28],  [29]  (see  also  [26]  for 
further  discussions  regarding  the  importance  of  such  problems). 

Before  using  our  estimation  or  inverse  techniques  in  problems  with  experimental  data, 
wc  carry  out  a  rather  careful  testing  of  the  methods  with  "synthetic"  data  on  numerous 
examples.  This  procedure  involves  a  series  of  tests  using  "data"  generated  (with  noise)  from 
a  system  with  known  (prechosen)  parameters  to  ascertain  the  ability  to  recover  the  parameters 
from  given  sets  of  "data".  For  detailed  explanations  of  this  procedure,  see  for  example  [3],  [5]. 
This  testing  is  also  combined  with  attempts  to  establish  convergence  and  stability  results  for  the 
methods.  For  problems  involving  systems  of  the  type  (3),  such  results  are  given  in  [1],  [10],  [II]. 
Wc  present  here  results  from  two  of  the  numerical  tests  we  performed. 


lie  1:  We  seek  to  estimate  D  =  D(u),  i.e.,  ?j,  a,  B  in  Figure  1,  in  the  system 

^  2u’-  +  F(t,x),  u(t,0)  -  u(t,l)  -  0,  u(0,x)  =  6x(l-x), 

at  ox  I  9xJ  2 


where  F  is  known  (computed  analytically  so  that  u(t,x)  =  6x(l-x)(l+t*)  is  a  solution 


corresponding  to  a  "true"  parameter  D*  with  to=.5,  {^=1.6,  a=3,  S=l).  "Data"  corresponding 
to  observations  at  (t^.Xj),  tpO,.5,l,Xj=.I,.2,...,.9.  were  used  for  the  inverse  procedure.  Results  for 
estimation  vkith  N=6  and  N=14  (cubic  splines  were  used  for  the  state  approximations  -  sec 
[10]  for  details)  along  with  the  initial  estimate  and  true  value  D*  are  depicted  in  Figure 
2. 


0.0  1.0  2.0 

Figure  2 

Example  2:  Wc  seek  to  estimate  q  and  the  nonlinearity  f  (we  do  not  make  any  a  priori 
parametrization  or  shape  assumption  on  f)  in 

=  q  +f(u)  +  F(t,x),  u(t,0)  =  u(t,l)  =  0,  u(0,x)  =  6x(l-x), 

9t 

where  again  F  is  computed  so  that  u(t,x)  =  6x{l-x)(l+t*)  is  a  solution  corresponding  to  true 
values  q*=  3.0,  f*(u)  =  2u^  -  u®.  Results  for  the  simultaneous  estimation  of  f  and  q  with 
state  approximation  N=6  (cubic  splines  with  mesh  h=  and  parameter  approximations 
(linear  splines)  for  f  with  mesh  size  h=.65  (sec  [II]  for  further  details)  are  presented  in 
Figure  3.  These  results  for  the  initial  estimate  f°,  the  converged  value  and  f* 

correspond  to  initial  guess  q°=1.0  and  converged  value  ^=2.9993. 


Figure  3 


SIZE  DEPENDENT  GROWTH  MODELS 

Wc  arc  currently  using  our  parameter  estimation  techniques  in  investigations  that  entail  size 
dependent  population  growth  models.  Data  from  experiments  with  mosquito  fish  populations  in 
rice  paddies  have  motivated  our  collaborative  efforts  with  L.  Botsford  (U.  California,  Davis). 
While  the  basic  problems  wc  are  considering  are  control  problems  for  the  mosquito 
fish/mosquito  populations,  a  substantial  effort  is  required  in  developing  the  underlying 
dynamics  (i.e.,  growth  models).  For  a  more  detailed  description  of  the  modeling  and  control 
problems,  wc  refer  to  [17],  [30]. 

A  simple  version  of  the  basic  modeling  problem  entails  estimation  of  q=(g,m,b)  in  the 
system 


dll 

at 


2_ 

d\ 


(gu) 


-mu 


XQ<x<Xj,t>0,  u(0,x)  =  fl^x),  g(xQ)u(t,Xo)=  J ’‘*b(t,Ou(t,Ocl?, 

*0 


(4) 


where  g=growth  rate,  m=mortality  rate,  and  b=fccundity  arc  in  general  dependent  on  time  t  and 
size  X,  with  x^,  Xj  the  minimum  and  maximum  observable  sizes,  respectively.  Data  for  the 
system  generally  consists  of  observations  that  yield  values  u(tj,x),  Xq«x«Xj,  so  that  a  distributed 
least-squares  criterion,  e.g.. 


J(q)  =  E  |u(t,.  ;q)-  u(tj,  )|  * 
i  L* 


is  appropriate. 


Since  a  convergence  theory  for  the  estimation  problem  has  not  appeared  elsewhere,  we 
sketch  one  approach  to  this  problem.  This  approach  is  the  analogue  to  that  given  in  [I],  [10], 
[II]  for  the  insect  dispersal  model  problems.  We  first  rewrite  (4)  in  variational  or  weak  form. 
We  seek  u(t)€H°(xQ,Xj)  satisfying  for  all  <i£H‘(Xg,Xj) 

+  <mu,4»  -  <gu,D4»  -  4>(xQ)/f(t,u)  =  0,  u(0, •)  =  ♦  (5) 

where  /?(t,ili)  =  f  *  b(t, Oi/^Odt,  D  =  ^,  and  <,>  is  the  usual  L*=H°  inner  product.  For  brevity, 

•**0  9* 

wc  assume  that  the  parameters  g  and  m  depend  only  on  size  x.  The  ideas  we  present  here  can 
be  readily  modified  to  treat  theoretically  and  computationally  the  more  general  case  where  g 
and  m  also  arc  time  dependent. 

Wc  assume  that  (g,m,b)  arc  to  be  chosen  from  function  spaces  CxMxffcH*(xQ,x  j)xH°(Xp, 
Xi)*H°([0,T]x[Xq,x  J)  containing  only  nonnegative  functions  g,m,b  where  the  functions  in  G  also 
satisfy  g(x  ^)=0  and  g(x^)>v^  for  some  positive  constant  Vj.  For  the  approximating  systems  (see 
[4],  [10])  we  assume  that  subspaccs  Z^CH‘(Xq,Xj)  are  chosen  and  let  u*^(t)£Z’^  denote  solutions  to 

<u[lil»  +  <mu^,t|»  -  <gu’^,D(lo'  -  iKxQ)/?(t,u^)=0,  for  all  u'^(0)=P^4>,  (6) 

where  is  the  orthogonal  projection  of  H°(Xq,Xj)  onto  Z^.  Wc  assume  that  the  subspaces 
satisfy: 

(Hj)  For  0£H‘  we  have  P^^-0  in  H‘  while  for  wc  have  in  H°. 

A  number  of  the  commonly  used  approximating  families  (piecewise  linear,  cubic  splines  [27]) 
satisfy  this  hypothesis.  The  ideas  here  can  be  slightly  modified  (see  the  remarks  in  [1])  to 
also  include  spectral  families  (such  as  Legendre  polynomials  -  sec  [20])  in  the  state 
approximation  schemes  for  which  the  convergence  theory  presented  here  is  valid. 

Wc  further  assume  that  <7,  A/,  B  arc  compact  in  H‘,  H°,  respectively,  and  that  compact 
approximation  families  C^,  b'^  for  G,  M,  B  respectively  have  been  chosen  satisfying: 

(Hj)  For  g'^(Xj)=0,  g’^(Xp)?Vj,  |Dg’^|*<jti,  and  C'^=i^(G)  where  for  each  geC,  i^(g)-g 

in  H*  with  the  convergence  uniform  in  geC; 

(Hj)  For  m^tA/^,  and  A/'^=i”(A/)  where  for  each  me  A/,  i5^(m)-m  in  H°  with  the 

convergence  uniform  in  meAf; 

(H^)  For  b*^£B*^,  lb^(t, ■  and  B^='i^(B)  where  for  each  beS,  i^(b)-b  in  H°  with  the 

convergence  uniform  in  btB. 


We  next  remark  that  to  give  a  convergence  theory  it  suffices  (see  [8])  to  argue  that 
uN(j.qM)_u({;q)  jn  N.M-®  whenever  [q'^)  is  an  arbitrary  sequence  with 

qMjj^Mx^/MxgM  gpjj  q^->qfGxAfxB.  Indeed,  it  suffices  to  give  these  arguments  in  the  form 
u*^(t;q^)-*u{t;q)  whenever  [q^)  is  arbitrary  in  CxA/xB  with  q^-*q  in  G^M*B  (in  the  H‘xH°xH° 
topology  in  this  case).  We  sketch  the  arguments;  let  (q^)  be  arbitrary  with  q^-q  in  Cx.t/xfl 
and  let  u'^^lq^),  u(q)  be  the  solutions  to  (6),  (5),  respectively,  corresponding  to  q^=(g*^,m^,b'^). 


q=(g,m,b)  respectively.  From  (Hj),  it  suffices  to  argue  that  z'^Ct)  h  u'^(t;q^)  -  P’^u(t;q)-0  in 
H°(Xq,Xj)  for  each  t  in  (0,T). 

Letting  =  /  *  b^(t,0<^<Od?.  we  have  from  (5)  and  (6)  that  for  all  iieZ^ 

*0 

<(uN-pNu)^  0>  =  <uN  .  +  (U-P’^'u)j,<» 

=  <(l-P’^)Uj,^>  +  <mu-m'^u^,^>  +  <g'^u'^  -  gu,D^>+  ^Xjj)A/{^' 
where  A/?^=^^(t,u*^)-V?(t,u).  With  as  defined  above  we  have  z’^(0)=0  and 
<z  ^,4»  =  <(l-P'^)Uj,^>  +  <mu-m^P^u,0>  -  <m^z^,4» 

+  <g^P^u-gu,Du>  +  <g'^z'^,D<»  +  ^Xq)A/?'^. 

Choosing  4>=2^  in  this  identity  we  find 


2  dt  '  ' 


=  <(I-P'^)Uj,z^>  +  <mu-m'^P’^u,z^>  +  <g^P'^u-gu,Dz'^>  +  z^(Xjj)A/?^. 
Recalling  that  g’^(Xj)  =  g(Xj)  =  0,  with  integration  by  parts  we  find 


<g'^z^Dz'^>  =  -  <Dg^z^z'^>  -  :^g^(X(j)  z^(Xo)^ 


and 


<g’^P^u-guX)z^>  =  -  <D(g^P^'u-gu),z^>  -  (g'^P^'u-gU](XQ)z'^(XQ). 
Hence,  we  have 


^  ^  <m^^z^'>  +  :!;■  <Dg^z^z^>  +  ^  g'^(Xo)z^(Xo)* 


=  <(I-P*^)Uj,z^>  +  <mu-m‘^P‘^u,z”>  -  <D(g'^P'^u-gu),z'^> 


NnN,.  ,N^ 


,NdN _ ■.  ,N^ 


(7) 


+  z'^'Cxq)  ((gu-g^P’^uKxg)  +  AR’^]. 

Using  the  bounds  from  (Hj),  (Hj)  we  find 

<m'^z^,z'^>  +  ^  <Dg'^z^,z^>  +  g^(XQ)z^(Xp)*  J  Hlz*^!^  +  vz'^(Xg)* 

for  positive  constants  (t,  v.  Standard  inequalities  imply  that  the  right  side  of  the  equality  (7)  is 
less  than  or  equal  to 


L_ 

4€ 


^|(I-P'^)Uj|  *  +  |mu-m^P^u|^  +  I  D(g^P'^u-gu)  I  ^  j 
+  3«  jz'^l*  +  €|z'^(Xo)|’  +  ^  |2|(gu-g'^p'^u)(Xo)|* 


where  €>0  is  arbitrarily  chosen.  From  the  bounds  of  (H^)  we  have  that  (/?''(t,z^)|«k|z^j  and 


defining  8'^(t)  =  /?^(t,P^u-u)  +  -  R(t,u)  we  may  conclude  that 

<  klz^l  +  |S^(t)j. 

Using  these  inequalities  in  (7)  we  obtain 

^  ^  |zN|2  +  |^4-3€-  |zN|2  +  (v-£)|zN(X(,)|*  <  h^d)  (8) 

where 

hN(t)  ||(I-P^)Uj|*  +  |mu-mNpNu|*  +  |D(g'‘P*'u  -  gu)|* 

+  2|(gu-g^P^u)(xQ)|*  +  4|6**(t)j*|. 


Choosing  €=v  and  using  z'^{0)=0,  we  may  apply  GronwalPs  lemma  to  conclude  that  it  suffices  to 
argue  that  h'^(t)-0  in  L‘(0,T)  to  obtain  the  desired  results.  Under  sufficient  smoothness  for  u(t), 
this  follows  readily  from  (H^)  -  (H^). 

While  we  arc  still  testing  our  methods  for  use  with  models  such  as  (4),  our  initial  findings 
are  quite  positive.  We  present  one  of  our  simple  test  examples. 

Example  .1:  We  seek  to  estimate  g  in  (4)  with  m=2,  b(x)=12x/i  -x,  4>(x)=»'2’t/j  - x,  Xq=0,  x,  =  1. 

The  solution  corresponding  to  g*(x)=2(l-x)  is  given  by  u(t,x)=eVr^i - x.  Cubic  splines  (Z^  of 
dimension  N+3)  were  used  for  the  state  approximations  while  linear  splines  (Q“  of  dimension 
M)  were  used  in  parameter  approximations  g^(x)  =  E  o.  I.  (x),  where  the  sum  is  from  i=0  to 
i=M-l  and  is  the  usual  "hat"  function  basis  element  (27]  with  support  in  ((i*l)/M,  (i+l)/M).  For 
M=4,  the  true  value  g*  corresponds  to  coefficients  o*,ot*)  =  (2.0,1.5,1.0,.5).  Using  data 

at  eleven  points  in  time  and  space  each  and  initial  guess  (1,1, 1,1),  we  obtained  the  converged 
values  (1.998,1.498,1.000,496)  for  a-cocfficients  in  the  representation  for  g^*'^  ,  i.c.,  with  N=32, 
M=4.  The  graphs  of  g*  and  g^*'*  are  not  distinguishable  using  ordinary  plotting  devices  and 
hence  we  do  not  present  them. 


COMPUTATIONAL  CONSIDERATIONS 


The  problems  on  which  we  have  focused  in  this  presentation  are  computationally 
intensive.  Even  simple  examples  such  as  those  presented  above  can  require  from  10^  to  10^ 
seconds  on  an  IBM  3081  and  we  are  now  using  the  ideas  discussed  here  in  research  problems  for 
which  use  of  such  a  sequential  machine  would  require  rather  prohibitive  computational 
expenditures.  The  necessary  software  packages  must  deal  with  reasonably  large  vector/tensor 
systems  and  involve  many  repetitive  routine  calculations.  Therefore,  a  substantial  part  of  our 
research  efforts  over  the  last  year  have  entailed  development  of  ideas,  algorithms,  and  software 
to  take  advantage  of  emerging  computer  architectures  involving  parallel  and  vector 
computational  capabilities.  Use  of  such  architectures  (in  our  research  programs  at  Brown 
University  and  ICASE,  we  are  currently  employing  a  CRAY  X-MP  -  a  widely  known  vector 
machine  -  and  a  STAR  ST-100  array  processor  with  parallel  features)  has  substantially  enhanced 


■  -  -  - 


our  efforts  to  investigate  some  of  the  research  questions  in  modeling  outlined  above.  For 
example,  results  for  problems  of  the  type  given  in  Example  1  typically  require  from  lO’  to  10^ 
seconds  on  the  IBM  3081,  but  when  the  algorithms  and  corresponding  software  are  modified  to 
take  advantage  of  the  arithmetic  speed  and  vector  capabilities  of  the  CRAY,  we  can  carry  out 
the  same  computational  runs  in  SO  to  200  seconds  on  a  CRAY  X-MP. 

Research  machines  such  as  the  CRAY  are,  through  NSF  and  other  research  sponsoring 
agencies,  becoming  readily  accessible  to  many  scientists  in  the  U.S.  We  recognize  that  the 
machines  we  are  presently  using  for  these  techniques  and  methods  are  not  widely  available  to 
the  world-wide  biological  research  community.  However,  we  firmly  believe  that  the  current 
revolution  in  computer  hardware  development  has  important  implications  for  the  community.  We 
note  that  many  of  the  high  speed,  parallel  and  vector  features  of  large,  expensive  research  array 
processors  such  as  the  FPS-164  and  STAR  ST-100  and  research  vector  machines  such  as  the 
CRAY  X-MP  (and  1-S)  and  CYBER  20S  are  rapidly  becoming  available  in  small,  relatively 
inexpensive  desk-top  configurations.  A  number  of  attached  array  processor  units  and  boards 
arc  now  available  for  use  with  personal  computers  such  as  the  IBM  PC  XT.  Recently,  new 
high  speed  chips  (e.g.,  INTEL  80386  -  a  32  bit,  4  MFLOPS  chip)  have  been  announced  and  will 
be  available  in  PC’s  in  1987.  We  believe  that  the  wide  availability  of  "desk-top  CRAY" 
capability  discussed  in  the  computer  science  community  is  only  several  years  away.  If  we  can 
successfully  develop  ideas  further  along  the  lines  of  those  discussed  above,  the  potential  for  a 
significant  impact  on  biological  modeling  and  research  is  enormous. 
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